#Amnesty Provision Analysis#

load("Amnesty Provision Main Data.Rdata")

#Figure 4 - Effect of Substitutes on Amnesty #####
#Continuous models  - w/ and w/o controls #
coefc<-matrix(nrow=1000, ncol=7, NA)
coefc_small<-matrix(nrow=1000,ncol=2, NA)
marginsc<-matrix(nrow=1000, ncol=2, NA)
databasec<-na.omit(data[,c("CID","Amn", "log_frag2")])
datafullc<-na.omit(data[,c("CID","Amn", "log_frag2", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                               "v2x_rule")])
popbase<-unique(databasec$CID)
popfull<-unique(datafullc$CID)
print("Starting Continuous Models")
set.seed(1015)
for(i in 1:1000){
  samplebase<-sample(popbase, replace = T, size = length(popbase))
  samplefull<-sample(popfull, replace = T, size = length(popfull))
  mapbase<-sapply(samplebase, function(x) which(databasec[,"CID"]==x))
  mapfull<-sapply(samplefull, function(x) which(datafullc[,"CID"]==x))
  boot.samplebase<-databasec[unlist(mapbase),]
  boot.samplefull<-datafullc[unlist(mapfull),]
  model1<-brglm(Amn~log_frag2,data=boot.samplebase, family="binomial")
  model2<-brglm(Amn~log_frag2+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull, family="binomial")
  coefc_small[i,]<-coef(model1)
  marginsc[i,1]<-summary(margins(model1))[1,2]
  coefc[i,]<-coef(model2)
  marginsc[i,2]<-summary(margins(model2))[2,2]
  #cat("\r", i)
}

#Binary models w/ and w/o controls#
coefb<-matrix(nrow=1000, ncol=7, NA)
coefb_small<-matrix(nrow=1000,ncol=2, NA)
marginsb<-matrix(nrow=1000, ncol=2, NA)
databaseb<-na.omit(data[,c("CID","Amn", "bin")])
datafullb<-na.omit(data[,c("CID","Amn", "bin", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                               "v2x_rule")])
popbase<-unique(databaseb$CID)
popfull<-unique(datafullb$CID)
set.seed(1015)
print("Starting Binary Models")
for(i in 1:1000){
  samplebase<-sample(popbase, replace = T, size = length(popbase))
  samplefull<-sample(popfull, replace = T, size = length(popfull))
  mapbase<-sapply(samplebase, function(x) which(databaseb[,"CID"]==x))
  mapfull<-sapply(samplefull, function(x) which(datafullb[,"CID"]==x))
  boot.samplebase<-databaseb[unlist(mapbase),]
  boot.samplefull<-datafullb[unlist(mapfull),]
  model1<-brglm(Amn~bin,data=boot.samplebase, family="binomial")
  model2<-brglm(Amn~bin+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull, family="binomial")
  coefb_small[i,]<-coef(model1)
  marginsb[i,1]<-summary(margins(model1))[1,2]
  coefb[i,]<-coef(model2)
  marginsb[i,2]<-summary(margins(model2))[1,2]
  #cat("\r", i)
}

#Pull together Figure 4 data (mfx)
#First, need the base mfx 

baseb<-brglm(Amn~bin, data=data)
fullb<-brglm(Amn~bin+Inc+v2x_clphy+v2x_cspart+v2x_polyarchy+
               v2x_rule, data=data)
mbaseb<-summary(margins(baseb))[1,2]
mfullb<-summary(margins(fullb))[1,2]

basec<-brglm(Amn~log_frag2, data=data)
fullc<-brglm(Amn~log_frag2+Inc+v2x_clphy+v2x_cspart+v2x_polyarchy+
               v2x_rule, data=data)
mbasec<-summary(margins(basec))[1,2]
mfullc<-summary(margins(fullc))[2,2]

plot<-rbind(mbaseb, mfullb, mbasec, mfullc)

ci<-rbind(t(apply(marginsb,2,sort)[c(25,975),]),
t(apply(marginsc,2,sort)[c(25,975),]))

plot<-cbind.data.frame(plot,ci)
colnames(plot)<-c("Estimate", "LB", "UB")
plot$Model<-c(2,2,1,1)
plot$Model<-factor(plot$Model, levels=c(2,1),
                      labels=c("Any Substitutes",
                               "Average Substitutes"))
plot$Controls<-c(1,2,1,2)
plot$Controls<-factor(plot$Controls, levels=c(1,2),
                           labels=c("No", "Yes"))
fig4<-ggplot(plot, aes(y=Estimate, x=reorder(Model, Estimate)))+
  coord_flip()+
  geom_point(aes(color=Controls), size=3, position=position_dodge(width=.5))+
  geom_linerange(aes(ymin=LB, ymax=UB, color=Controls), position=position_dodge(width=.5))+
  geom_hline(yintercept=0)+
  xlab("")+
  scale_color_manual(values=viridis_pal()(3)[c(1,2)])+
  ylab("Effect of Substitutes on Amnesty Provision")+
  theme_bw()
fig4
pdf("Fig4.pdf", width=6, height=2)
print(fig4)
dev.off()


#Appendix Table 4 - Effect of Substitutes on Amnesty ######
#Model 1 - Continuous, No Controls
m1<-cbind(coef(basec), 
      t(apply(coefc_small,2,sort)[c(25,975),]))
colnames(m1)<-c("Estimate", "Lower Bound", "Upper Bound")
names(plot)<-c("Estimate", "Lower Bound", "Upper Bound")
m1<-rbind(m1, plot[3,1:3])
m1<-data.frame(m1)
m1<-round(m1,3)
m1$CI<-paste("[", m1$Lower.Bound, " --- ", m1$Upper.Bound, "]", sep="")
m1<-as.vector(t(m1[,c("Estimate", "CI")]))
m1[length(m1)+1]<-nobs(basec)
m1[length(m1)+1]<-round(logLik(basec), digits=3)
m1


m2<-cbind(coef(fullc), 
          t(apply(coefc,2,sort)[c(25,975),]))
colnames(m2)<-c("Estimate", "Lower Bound", "Upper Bound")
m2<-rbind(m2, plot[4,1:3])
m2<-data.frame(m2)
m2<-round(m2,3)
m2$CI<-paste("[", m2$Lower.Bound, " --- ", m2$Upper.Bound, "]", sep="")
m2<-as.vector(t(m2[,c("Estimate", "CI")]))
m2[length(m2)+1]<-nobs(fullc)
m2[length(m2)+1]<-round(logLik(fullc), digits=3)
m2


m3<-cbind(coef(baseb), 
          t(apply(coefb_small,2,sort)[c(25,975),]))
colnames(m3)<-c("Estimate", "Lower Bound", "Upper Bound")
m3<-rbind(m3, plot[1,1:3])
m3<-data.frame(m3)
m3<-round(m3,3)
m3$CI<-paste("[", m3$Lower.Bound, " --- ", m3$Upper.Bound, "]", sep="")
m3<-as.vector(t(m3[,c("Estimate", "CI")]))
m3[length(m3)+1]<-nobs(baseb)
m3[length(m3)+1]<-round(logLik(baseb), digits=3)
m3

m4<-cbind(coef(fullb), 
          t(apply(coefb,2,sort)[c(25,975),]))
colnames(m4)<-c("Estimate", "Lower Bound", "Upper Bound")
m4<-rbind(m4, plot[2,1:3])
m4<-data.frame(m4)
m4<-round(m4,3)
m4$CI<-paste("[", m4$Lower.Bound, " --- ", m4$Upper.Bound, "]", sep="")
m4<-as.vector(t(m4[,c("Estimate", "CI")]))
m4[length(m4)+1]<-nobs(fullb)
m4[length(m4)+1]<-round(logLik(fullb), digits=3)
m4

#Get everything to same length in right order
m1<-as.matrix(m1)
#Move intercept to bottom of coef
m1<-as.matrix(m1[c(3:4,1:2,5:8),])
#M1 needs 12 blanks in 3:15
m1<-rbind(cbind(m1[1:2]), cbind(rep(NA,12)),
      cbind(m1[3:nrow(m1)]))
m2<-as.matrix(m2)
m2<-as.matrix(m2[c(3:14,1:2,15:length(m2)),])
#M2 needs 2 blanks in 3:4
m2<-rbind(cbind(m2[1:2]),cbind(rep(NA, 2)),
            cbind(m2[3:nrow(m2)]))
m3<-as.matrix(m3)
#Move intercept
m3<-as.matrix(m3[c(3:4,1:2,5:8),])
#M3 needs two blanks at the top,
#then 10 in the 5:14 position
m3<-rbind(cbind(rep(NA,2)),
      cbind(m3[1:2]),
      cbind(rep(NA,10)),
      cbind(m3[3:nrow(m3)]))
#Move intercept
m4<-as.matrix(m4)
m4<-as.matrix(m4[c(3:14,1:2,15:length(m4)),])
#M4 needs two blanks on top and that's it
m4<-rbind(cbind(rep(NA,2)),
      cbind(m4))
all<-cbind.data.frame(m1,m2,m3,m4)
rownames(all)[c(seq(1,17,2),19:20)]<-c("Average Substitutes",
                "Any Substitutes",
                "Incompatibility",
                "Political Violence",
                "Civil Society Strength",
                "Electoral Democracy",
                "Rule of Law",
                "Intercept",
                "Marginal Effect of Substitutes",
                "N",
                "Log-Likelihood")
#names(all)<-c("")
addtorow <- list()
addtorow$pos <- list(0)
addtorow$command <- "& \\multicolumn{2}{c}{Average Substitutes} & \\multicolumn{2}{c}{Any Substitutes}\\\\\n" 
all<-print(xtable(all, digits=3, align=c("l","c","c","c","c")), hline.after = c(16,18), include.colnames = F,
           add.to.row = addtorow)
write(all, file="Table 4.tex")


#Appendix Figure 6 - Alternative Measures of Substitutes #####
#Get point estimates from all 8 model specifications
basec1_no<-brglm(Amn~log(fragmentation2_noinfight+1), data=data)
basec2_no<-brglm(Amn~log(fragmentation2_noinfight+1) +as.factor(Inc)+
                   v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule , data=data)
basec1_forge<-brglm(Amn~log(fragmentation2_forge+1), data=data)
basec2_forge<-brglm(Amn~log(fragmentation2_forge+1) +as.factor(Inc)+
                       v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule , data=data)

baseb1_no<-brglm(Amn~bin_noinfight, data=data)
baseb2_no<-brglm(Amn~bin_noinfight +as.factor(Inc)+
                   v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule , data=data)
baseb1_forge<-brglm(Amn~bin_forge, data=data)
baseb2_forge<-brglm(Amn~bin_forge +as.factor(Inc)+
                       v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule , data=data)

#Begin cluster bootstrap

marginsb<-matrix(nrow=1000, ncol=4, NA)

databaseb_no<-na.omit(data[,c("CID","Amn", "bin_noinfight")])
datafullb_no<-na.omit(data[,c("CID","Amn", "bin_noinfight", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                                  "v2x_rule")])
databaseb_s<-na.omit(data[,c("CID","Amn", "bin_forge")])
datafullb_s<-na.omit(data[,c("CID","Amn", "bin_forge", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                                 "v2x_rule")])

popbase_no<-unique(databaseb_no$CID)
popfull_no<-unique(datafullb_no$CID)
popbase_s<-unique(databaseb_s$CID)
popfull_s<-unique(datafullb_s$CID)

set.seed(1015)
#All Binary Measures#
for(i in 1:1000){
  samplebase_no<-sample(popbase_no, replace = T, size = length(popbase_no))
  samplefull_no<-sample(popfull_no, replace = T, size = length(popfull_no))
  samplebase_s<-sample(popbase_s, replace = T, size = length(popbase_s))
  samplefull_s<-sample(popfull_s, replace = T, size = length(popfull_s))
  
  mapbase_no<-sapply(samplebase_no, function(x) which(databaseb_no[,"CID"]==x))
  mapfull_no<-sapply(samplefull_no, function(x) which(datafullb_no[,"CID"]==x))
  mapbase_s<-sapply(samplebase_s, function(x) which(databaseb_s[,"CID"]==x))
  mapfull_s<-sapply(samplefull_s, function(x) which(datafullb_s[,"CID"]==x))
  
  boot.samplebase_no<-databaseb_no[unlist(mapbase_no),]
  boot.samplefull_no<-datafullb_no[unlist(mapfull_no),]
  boot.samplebase_s<-databaseb_s[unlist(mapbase_s),]
  boot.samplefull_s<-datafullb_s[unlist(mapfull_s),]
  
  model1<-brglm(Amn~bin_noinfight,data=boot.samplebase_no, family="binomial")
  model2<-brglm(Amn~bin_noinfight+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull_no, family="binomial")
  marginsb[i,1]<-summary(margins(model1))[1,2]
  marginsb[i,2]<-summary(margins(model2))[1,2]
  
  model3<-brglm(Amn~bin_forge,data=boot.samplebase_s, family="binomial")
  model4<-brglm(Amn~bin_forge+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull_s, family="binomial")
  marginsb[i,3]<-summary(margins(model3))[1,2]
  marginsb[i,4]<-summary(margins(model4))[1,2]
  #cat("\r", i)
}

#All Continuous Measures #

marginsc<-matrix(nrow=1000, ncol=4, NA)
databaseb_no<-na.omit(data[,c("CID","Amn", "fragmentation2_noinfight")])
datafullb_no<-na.omit(data[,c("CID","Amn", "fragmentation2_noinfight", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                                  "v2x_rule")])
databaseb_s<-na.omit(data[,c("CID","Amn", "fragmentation2_forge")])
datafullb_s<-na.omit(data[,c("CID","Amn", "fragmentation2_forge", "Inc","v2x_cspart","v2x_polyarchy","v2x_clphy",
                                 "v2x_rule")])
set.seed(1015)
for(i in 1:1000){
  
  samplebase_no<-sample(popbase_no, replace = T, size = length(popbase_no))
  samplefull_no<-sample(popfull_no, replace = T, size = length(popfull_no))
  samplebase_s<-sample(popbase_s, replace = T, size = length(popbase_s))
  samplefull_s<-sample(popfull_s, replace = T, size = length(popfull_s))
  
  mapbase_no<-sapply(samplebase_no, function(x) which(databaseb_no[,"CID"]==x))
  mapfull_no<-sapply(samplefull_no, function(x) which(datafullb_no[,"CID"]==x))
  mapbase_s<-sapply(samplebase_s, function(x) which(databaseb_s[,"CID"]==x))
  mapfull_s<-sapply(samplefull_s, function(x) which(datafullb_s[,"CID"]==x))
  
  boot.samplebase_no<-databaseb_no[unlist(mapbase_no),]
  boot.samplefull_no<-datafullb_no[unlist(mapfull_no),]
  boot.samplebase_s<-databaseb_s[unlist(mapbase_s),]
  boot.samplefull_s<-datafullb_s[unlist(mapfull_s),]
  
  model1<-brglm(Amn~log(fragmentation2_noinfight+1),data=boot.samplebase_no, family="binomial")
  model2<-brglm(Amn~log(fragmentation2_noinfight+1)+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull_no, family="binomial")
  marginsc[i,1]<-summary(margins(model1))[1,2]
  marginsc[i,2]<-summary(margins(model2))[1,2]
  
  model3<-brglm(Amn~log(fragmentation2_forge+1),data=boot.samplebase_s, family="binomial")
  model4<-brglm(Amn~log(fragmentation2_forge+1)+as.factor(Inc)+
                  v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=boot.samplefull_s, family="binomial")
  marginsc[i,3]<-summary(margins(model3))[1,2]
  marginsc[i,4]<-summary(margins(model4))[1,2]
  #cat("\r", i)
}

plot<-rbind.data.frame(t(apply(marginsc,2,sort)[c(25,975),]),
                       t(apply(marginsb,2,sort)[c(25,975),]))
point<-rbind(summary(margins(basec1_no))[1,2],
             summary(margins(basec2_no))[1,2],
             summary(margins(basec1_forge))[1,2],
             summary(margins(basec2_forge))[1,2],
             summary(margins(baseb1_no))[1,2],
             summary(margins(baseb2_no))[1,2],
             summary(margins(baseb1_forge))[1,2],
             summary(margins(baseb2_forge))[1,2])

plot$point<-point
plot$measure<-c(rep("No Infighting (Continuous)",2),
                rep("FORGE (Continuous)", 2),
                rep("No Infighting (Binary)",2),
                rep("FORGE (Binary)", 2))
plot$Controls<-c(rep(c("No", "Yes"),4))

plot$ind<-rep(c(1.5,2),4)
p<-ggplot(plot, aes(y=point, x=ind))+
  coord_flip()+
  geom_point(aes(color=Controls), size=3, position=position_dodge(width=.5))+
  geom_linerange(aes(ymin=V1, ymax=V2, color=Controls), position=position_dodge(width=.5))+
  geom_hline(yintercept=0)+
  scale_x_continuous("", limits=c(0,4), breaks=NULL)+
  ylab("Effect of Substitutes on Amnesty Provision")+
  scale_color_manual(values=viridis_pal()(3)[c(1,2)])+
  facet_wrap(~measure)+
  theme_bw()
p
pdf("FigA6.pdf", height = 4, width = 5)
print(p)
dev.off()

#Appendix Figure 7 - Group Level Controls #####
basec<-brglm(Amn~log_frag2 +as.factor(Inc)+
               v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule+
               groupcont+groupcentany+groupsize +groupresource, data=data)
baseb<-brglm(Amn~bin +as.factor(Inc)+
               v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule+
               groupcont+groupcentany+groupsize +groupresource, data=data)

databaseb<-na.omit(data[,c("CID","Amn", "bin", "log_frag2", "Inc",
                               "v2x_cspart","v2x_polyarchy", "v2x_clphy",
                               "v2x_rule", "groupcont","groupcentany","groupsize","groupresource")])

popbaseb<-unique(databaseb$CID)

margins<-matrix(nrow=1000, ncol=2, NA)
set.seed(1015)
i<-0
n<-1000
while(length(na.omit(margins[,2]))<n){
  i<-i+1
  holder<-matrix(NA, ncol=2)
  samplebaseb<-sample(popbaseb, replace = T, size = length(popbaseb))
  
  mapbaseb<-sapply(samplebaseb, function(x) which(databaseb[,"CID"]==x))
  
  boot.samplebaseb<-databaseb[unlist(mapbaseb),]
  #Kick out if non-varying on any X
  if(min(apply(boot.samplebaseb, 2, var))==0) 
    next
  #Kick out if singular vcov mat
  chol<-NA
  try(chol<-chol.default(cov(boot.samplebaseb)))
  if(any(is.na(chol))) next
  try(model1<-brglm(Amn~bin+as.factor(Inc)+
                      v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule+
                      groupsize+groupcont+groupcentany+groupresource, 
                    data=boot.samplebaseb, family="binomial"))
  if(model1$converged==FALSE) next
  try(model2<-brglm(Amn~log_frag2+as.factor(Inc)+
                      v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule+
                      groupsize+groupcont+groupcentany+groupresource, 
                    data=boot.samplebaseb, family="binomial"))
  if(model2$converged==FALSE) next
  holder[,1]<-summary(margins(model1))[1,2]
  holder[,2]<-summary(margins(model2))[6,2]
  margins<-rbind(margins,holder)
  #cat("\r", i)
}


plot<-data.frame(t(apply(margins,2,sort)[c(25,975),]))
plot$point<-NA
plot[1,]$point<-summary(margins(baseb))[1,2]
plot[2,]$point<-summary(margins(basec))[6,2]
plot$measure<-c("Binary", "Continuous")
colnames(plot)[1:2]<-c("lb","ub")

fig7<-ggplot(plot, aes(y=point, x=measure))+
  coord_flip()+
  geom_point(size=3, position=position_dodge(width=.5))+
  geom_linerange(aes(ymin=lb, ymax=ub), position=position_dodge(width=.5))+
  geom_hline(yintercept=0)+
  ylab("Effect of Substitutes on Amnesty Provision")+
  xlab("Measure")+
  theme_bw()
fig7
pdf("Fig7.pdf", height = 4, width = 4)
print(fig7)
dev.off() 
#Appendix Figures 11 and 12 - Including One Control #####
varlist<-c("as.factor(Inc)","v2x_cspart", "v2x_polyarchy",
           "v2x_clphy", "v2x_rule")
plotc<-NA
plotb<-NA

#Double loop - inner runs the bootstrap 1000 times
#Outer iterates over the control list to include one
#This takes a long time
set.seed(1015)
for (j in 1:5){
  formulac<-as.formula(paste0("Amn~log_frag2+", paste0(varlist[j])))
  formulab<-as.formula(paste0("Amn~bin+", paste0(varlist[j])))
  
  basec<-brglm(formulac, data=data)
  baseb<-brglm(formulab, data=data)
  
  coefc<-matrix(NA, nrow=1, ncol=1000)
  coefb<-matrix(NA, nrow=1, ncol=1000)
  for(i in 1:1000){
    sample<-sample(unique(data$CID), 
                   replace = T, size = length(unique(data$CID)))
    map<-sapply(sample, function(x) which(data[,"CID"]==x))
    boot.sample<-data[unlist(map),]
    bootc<-brglm(formulac, data=boot.sample)
    bootb<-brglm(formulab, data=boot.sample)
    
    coefc[i]<-summary(margins(bootc))[which(summary(margins(basec))$factor==
                                              "log_frag2"),2]
    coefb[i]<-summary(margins(bootb))[which(summary(margins(baseb))$factor==
                                              "bin"),2]
    #cat("\r", paste(j, i, sep=":"))
  }
  plotc<-rbind(plotc,cbind.data.frame(varlist[j],summary(margins(basec))[which(summary(margins(basec))$factor==
                                                                                 "log_frag2"),2],
                                      sort(coefc, decreasing=F)[25], sort(coefc, decreasing=F)[975]))
  
  plotb<-rbind(plotb,cbind.data.frame(varlist[j],summary(margins(baseb))[which(summary(margins(baseb))$factor==
                                                                                 "bin"),2],
                                      sort(coefb, decreasing=F)[25], sort(coefb, decreasing=F)[975]))
  
}  


plotc<-na.omit(plotc)
plotb<-na.omit(plotb)

colnames(plotc)<-c("Variable", "Point", "LB", "UB")
colnames(plotb)<-c("Variable", "Point", "LB", "UB")

plotc$Label<-seq(1,5,1)
plotc$Label<-factor(plotc$Label,
                    labels=c("Incompatibility",
                             "Civil Society Strength",
                             "Electoral Democracy",
                             "Political Violence",
                             "Rule of Law"))

fig11<-ggplot(plotc, aes(y=Point, x=Label))+
        coord_flip()+
        geom_point(size=3, position=position_dodge(width=.5))+
        geom_linerange(aes(ymin=LB, ymax=UB), position=position_dodge(width=.5))+
        geom_hline(yintercept=0)+
        xlab("Included Control")+
        ylab("Effect of Substitutes on Amnesty Provision")+
        theme_bw()
fig11
pdf("Fig11.pdf",width=6, height=4)
print(fig11)
dev.off()

plotb$Label<-seq(1,5,1)
plotb$Label<-factor(plotb$Label,
                    labels=c("Incompatibility",
                             "Civil Society Strength",
                             "Electoral Democracy",
                             "Political Violence",
                             "Rule of Law"))
fig12<-ggplot(plotb, aes(y=Point, x=Label))+
        coord_flip()+
        geom_point(size=3, position=position_dodge(width=.5))+
        geom_linerange(aes(ymin=LB, ymax=UB), position=position_dodge(width=.5))+
        geom_hline(yintercept=0)+
        xlab("Included Control")+
        ylab("Effect of Substitutes on Amnesty Provision")+
        theme_bw()
fig12
pdf("Fig12.pdf",width=6, height=4)
print(fig12)
dev.off()

#Appendix Figures 13 and 14 - Bayesian Model Averaging#####
bma_dat<-na.omit(data[,c("Amn", "log_frag2", "bin", "Inc", "v2x_cspart", "v2x_polyarchy", "v2x_clphy", "v2x_rule")])

colnames(bma_dat)<-c("Amnesty Provision",
                 "Average Substitutes",
                 "Any Substitutes",
                 "Incompatibility",
                 "Civil Society Strength",
                 "Electoral Democracy",
                 "Political Violence",
                 "Rule of Law")

#Figure 13 - Continuous Measure of Substitutes
bma_datc<-bma_dat[,(names(bma_dat)!="Any Substitutes")]
set.seed(1015)
bas<-bas.glm(`Amnesty Provision`~., data=bma_datc, family="binomial"(link="logit"), betaprior = g.prior(3))

plot(coef(bas), subset=2, ask=F)

pdf("Fig13.pdf", width = 6, height = 3)
print(plot(coef(bas), subset=2, ask=F))
dev.off()

#Figure 14- Binary Measure
bas_datb<-bma_dat[,(names(bma_dat)!="Average Substitutes")]
set.seed(1015)
bas2<-bas.glm(`Amnesty Provision`~., data=bas_datb, family="binomial"(link="logit"), betaprior = g.prior(3))

plot(coef(bas2), subset=2, ask=F)

pdf("Fig14.pdf", width = 6, height = 3)
print(plot(coef(bas2), subset=2, ask=F))
dev.off()

rm(list=setdiff(ls(), c("cluster_bootstrap", "int_plot")))